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ABSTRACT 

We revisit the idea that density-wave wakes of planets drive accretion in protostel- 
lar disks. The effects of many small planets can be represented as a viscosity if the 
wakes damp locally, but the viscosity is proportional to the damping length. Damping 
occurs mainly by shocks even for earth-mass planets. The excitation of the wake fol- 
lows from standard linear theory including the torque cutoff. We use this as input to 
an approximate but quantitative nonlinear theory based on Burger's equation for the 
subsequent propagation and shock. Shock damping is indeed local but weakly so. If all 
metals in a minimum-mass solar nebula are invested in planets of a few earth masses 
each, dimensionless viscosities [alpha] of order dex(-4) to dex(-3) result. We compare 
this with observational constraints. Such small planets would have escaped detection 
in radial-velocity surveys and could be ubiquitous. If so, then the similarity of the 
observed lifetime of T Tauri disks to the theoretical timescale for assembling a rocky 
planet may be fate rather than coincidence. 

Subject headings: planets and satellites: general — solar system: formation — (stars:) 
planetary systems 



1. Introduction 



The disks of young low- mass protostars ( "protostellar disks") are observed to accrete at rates 
~ 1CP 8 Mq yr _1 (Hartmann et al. 1998). A comparable average accretion rate is implied by the 
disk masses inferred from dust emission (~ 10~ 2±1 Mq, Osterloh & Beckwith 1995) when combined 
with the maximum ages of T Tauri stars that show evidence for disks (~ 10 6-7 yr, Strom et al. 
1989), or with radiochemical estimates for the lifetime of the protostellar nebula (Podosek & Cassen 
1994). Indeed, the evidence that young protostars are surrounded by viscous accretion disks has 
only strengthened in the quarter century since Lynden-Bell & Pringle (1974) proposed it. 

While the effective viscosity of the hotter accretion disks surrounding degenerate stars and 
black holes is probably due to magnetohydrodynamic (MHD) turbulence (cf. Balbus & Hawley 
1998), protoplanetary disks may be too cold, too weakly ionized, and too resistive to sustain MHD 
turbulence (Blaes & Balbus 1994; Jin 1996; Hawley & Stone 1998) except perhaps at small radii, 
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in surface layers (Gammie 1996), or during FU Orionis-type outbursts (Hartmann et al. 1993). In 
view of several physical uncertainties in the ionization balance (cf. Gammie 1996), and the possibly 
subtle nature of the conductivity itself (Wardle 1999), this conclusion is only tentative. But it is 
probably fair to say that the best reason for supposing that the viscosity of protoplanetary disks is 
magnetic remains the lack of a theoretically attractive alternative. Convection (Lin & Papaloizou 
1980) and nonlinear shear instability (Zahn 1991; Dubrulle 1992) have been proposed, but these 
mechanisms appear to fail because of the very strongly stabilizing influence of a positive radial 
gradient in specific angular momentum (Ryu & Goodman 1992; Stone & Balbus 1996; Balbus et 
al. 1996). 

Larson (1989) suggested that density waves might provide the efffective viscosity mechanism 
in protostellar disks. As he made clear, wave transport is rather different from viscous transport. 
In the standard theory of thin accretion diks, viscosity (y) plays two roles. Firstly, it provides an 
outward angular-momentum flux (Fj), a torque exerted by the part of the disk interior to a given 
radius on the part exterior. Secondly, it heats the disk by dissipating the differential rotation. In 
steady state, Fj balances the inward flux of angular momentum carried by the accreting matter, 
while the dissipation balances the loss of its orbital energy. Even if the disk is not steady, provided 
only that it is truly viscous, then 

n dQ, 

Fj = -u x 27rEr 3 V 
dr 

where S is the surface density and Q the angular velocity. Also, the dissipation rate per unit area 
is related to Fj by 

dQ 

E = -Fj-. (1) 

Any dissipative mechanism such that Fj cx d^l/dr and eq. (1) hold can be described by an effective 
viscosity u e g; one expects this to be the case for dissipative mechanisms that are sufficiently local. 

Since density waves are not local, one may deduce a different value of v e g from Fj than from E. 
In particular, inviscid linear theory (Goldreich & Tremaine 1980, henceforth GT80) predicts that a 
planet embedded in a gas disk excites density waves at the inner and outer Lindblad resonances of 
its (circular) orbit; these waves carry Fj > but E = by assumption, so that = but vj > 0. 
Actually, the theory implicitly assumes some damping at large distance from the resonances, else 
the waves would reflect from the edges of the disk and produce zero net flux. If one is concerned 
with the torques on the disk only, then the damping length is unimportant as long as it is large. 
But without dissipation, the waves cause no secular changes to the orbits of the gas even at the 
Lindblad resonances (Goldreich & Nicholson 1989). 

By demanding that eq. (1) be satisfied, Larson (1990) found that steady spiral shocks produce 

a ~ 0.013 [{c/Vf + 0.08(c/F) 2 ] 1/2 (2) 



independently of the source of excitation of the waves. Here c is the sound speed and V = rQ the 
orbital velocity. For c/V = 0.03, perhaps characteristic of the later phases of the protosolar nebula 
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at 1 AU, this gives a 10~ 4 . Larson's result depends upon a number of additional assumptions, 
some inspired by earlier work by Spruit (1987) on self-similar shocks: the waves were taken to form 
two spiral arms separated by 180° in azimuth; and because they would be very tightly wrapped since 
c<y, Larson (1990) approximated the wave profile by analyzing a train of perfectly axisymmetric 
nonlinear waves. Most importantly, at least by contrast with our own work reported below, no 
source of excitation for the waves exists within the range of radii where eq. (2) applies. Larson 
argued, however, that a bisymmetric spiral shock of the required strength could be launched from 
the 2 : 1 resonance of a jovian-mass planet, and he deduced an accretion timescale ~ 10 6-7 yr in 
accord with observational requirements. 

Jupiter-mass planets have been discovered during the past few years by exquisitely sensitive 
radial- velocity surveys (Marcy et al. 2000). Most of these planets have orbital periods shorter than 
that of Mercury. They are detected around 4 — 5% of the stars surveyed, mostly G, F, and K 
dwarfs. There are indications that stars of higher-than-solar metallicity are preferred (Gonzalez 
1997; Santos et al. 2000; Gonzalez et al. 2000). Because of observational selection effects, it is 
probably too early to derive the incidence of giant planets in long-period orbits (~ 10 yr). If such 
planets are as rare as those on short-period orbits, however, then we must seek another explanation 
for the observed accretion rates and inferred lifetimes of most protostellar disks. Nevertheless, it 
remains an intriguing notion that planet formation is not just constrained by disk lifetimes, but 
may actually determine them. 

Our paper therefore inquires whether accretion could be driven by planets of roughly terrestrial 
mass, which radial-velocity surveys would not have detected. Other things being equal, the torque 
exerted on the disk by a planet is proportional to the square of the planet's mass (M p ), so the 
influence of M p ~ 1O~ 5 M might be thought negligible compared with that of M p ~ 10 3 M Q . 
Just because its local effects are so strong, however, the larger planet will almost certainly open 
a gap in the disk, whereas the smaller one will not; since the torque per unit radius (smoothed 
over resonances) decreases rapidly with distance from the planet (GT80), the shorter range of its 
interaction may partly compensate for the terrestrial planet's smaller mass. At solar abundance 
there could be tens of earth-mass planets in the disk, each too small to have captured a significantly 
massive atmosphere. We assume that the planets are distributed through the disk in proportion to 
the local surface density of the gas, neglecting orbital migration and radial changes in the timescale 
of planet formation. Unlike Larson, therefore, we have a distributed source of excitation for density 
waves, so that eq. (2) need not apply (although in fact we obtain a comparable value of a in the 
minimum-mass solar nebula). 

As discussed, no sensible effective viscosity can be derived for density waves without considering 
their damping. §2 shows for a generic damping mechanism in a thin disk that if the waves are 
absorbed over a distance small compared to the disk radius, but possibly large compared to its 
thickness, and if the local damping rate depends mainly on distance from the planet, then eq. (1) 
is satisfied so that a consistent value of v e & or a can in fact be defined. However, the viscosity is 
directly proportional to the damping length. After briefly considering other mechanisms in §3, we 
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conclude that shock damping is probably most important even at M p ~ M®. In §4, we calculate 
the excitation of the wake from linear theory, which is valid for sufficiently small M p [eq. (31]. In §5 
and the Appendix, we show that the subsequent nonlinear propagation of the wake away from the 
planet can be reduced to Burger's equation using some transformations and some approximations 
that should be accurate for sufficiently low-mass planets in sufficiently thin disks. We describe 
quantitatively the formation of the shock and its subsequent decay with increasing distance from 
the planet. Most importantly, we obtain formulae for the shock-damping length, and hence for a e g, 
as functions of M p and local disk parameters. §6 summarizes our main results and discusses their 
implications for disks resembling a standard minimum-mass solar nebula. 



2. Effective a-parameter. 

Let f^p (M p ,r p ) be the angular- momentum flux excited by a planet of small mass M p and 
orbital radius r p , embedded in a disk of local surface density S, sound speed c, and angular velocity 
£l(r) ~ (GM^/r 3 ) 1 / 2 , where M* is the mass of the central star. According to GT80, 

ff\M p ,r p ) = {GM p f^ x \^nl ax {Q) [2*0(2/3) + *i(2/3)] 2 } . (3) 

The dimensionless quantity [i is defined in such a way that fj,Clr/c 3> 1 is roughly the largest az- 
imuthal harmonic contributing to the flux (the "torque cutoff" ) ; it is a function of the gravitational 
stability parameter Q = ctl/irGY, (Toomre 1964). Selfgravity is probably very weak in the later 
evolutionary phases of protostellar disks, so it is appropriate to take Q — > oo. In this limit, it 
follows from GT80's results that the quantity in curly brackets ~ 0.93. As acknowledged by GT80 
and elaborated by later authors, eq. (3) is valid only to leading order in c/Clr, the ratio of disk 
thickness to radius. At the next order, there is a slight imbalance between the interior (r < r p ) 
and exterior (r > r p ) fluxes; the difference between the two torques represents a net torque on the 
planet, whose orbit therefore "migrates," probably towards smaller r p (Ward 1997, and references 
therein). The effective viscosity we calculate depends on the sum of the interior and exterior fluxes, 
so we neglect the difference between them. 

The flux rises quickly to a constant value (3) over a radial distance from the planet of order 
I = c/\rd£l/dr\~ l , the length over which the rotation velocity changes by Mach 1. The strongest 
Lindblad resonances giving rise to the wake lie within this distance. Absent damping, the angular 
momentum flux is radially constant at |r — r p \ 3> I. We can characterize the strength of damping 
by the reciprocal of the damping length Id at which fj falls to ~ e _1 of eq. (3). 

The larger the damping length, the greater the dissipation of orbital energy by the density- 
wave wake, for the following reason. If a quantity A J of angular momentum is added to the disk 
in the form of waves having angular pattern speed Q, p (in steady state, the wake has the pattern 
speed of the planet), the work required to create the waves is Q p AJ. When these waves damp at 
radius r, the work done on the mean flow there is Q(r) A J. The difference [f2 p — Q(r)]AJ represents 
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the energy available to be dissipated, assuming that this is positive. It can be shown that AJ 
is positive (negative) for waves excited exterior (interior) to the planet in the approximation that 
Sip = £l(r p ). From these considerations, it follows that the total mechanical energy dissipated by 
the wake of the one planet is 

oo 

e(M p ,r p ) = J [n(r)-n p ]^-dr. (4) 
o 

Clearly e oc Id, asssuming Id <C r p . Suppose that the background properties of the disk are 
approximately constant over radial distances of order Id- Write the angular momentum flux in terms 
of its ideal asymptotic value (3) and a dimensionless distribution function, fj(r) = ff^ ■ 4>(r — r p ). 
Our notation presumes a damping mechanism such that 4> depends much more rapidly on r — r p 
than on (r + r p )/2, the latter dependence being significant only on the scale of r itself. The integral 
(4) can now be approximated by 

oo 

e(M p ,r p ) = f?\M p ,r p )(fy J xfjx, (5) 



— oo 



where x = r — r p . The lower limit of integration has been extended to — oo for convenience since 
4>{x) ps unless \x\ <C r p . 

Now suppose that there are many planets of the same mass distributed through the disk in 
numbers N(r) per unit radius, and that rN(r) 3> 1. Then the total flux at radius r is (on average) 

oo oo 

Fj(r) = J N(r p )ff\M p ,r p )<P(r-r p )dr p « N{r)ff{M p , r) j <f>(x)dx, (6) 

-oo 

assuming that N(r) and ff\M p ,r) vary only slowly with radius. Similarly, the dissipation rate 
per unit radius at r is 

oo 

E(r) = j[n(r)-n(r p )]N(r p )ff\M p ,r p )^(r-r p )dr p 
o 

oo 

H £*(rWf(M„r) / = -§FAr), (7) 

— oo 

where the last step has used integration by parts. So, as promised in §1, equation (1) is satisfied 
by density-wave wakes under the assumptions of many planets and local damping. 

In the standard theory of viscous accretion disks (Pringle 1982), the viscous energy dissipation 
rate per unit radius is related to the Shakura & Sunyaev (1973) viscosity parameter a by (assuming 
n oc r" 3 / 2 ) 

E visc = a—r^c 2 . (8) 
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Let the planets have average mass (M p ) and let them account for a fraction Z v <C 1 of the total 
nebular mass; we expect that Z v is comparable to the metallicity of the gas if planet-formation is 
efficient but the planets are too small to have captured massive gaseous envelopes (Mizuno et al. 
1978; Podolak et al. 1993; Ikoma et al. 2000) Then the average number of planets per unit radius 

2vrrZ p E 

(Mp) ' ( j 

The total dissipation rate per unit radius due to planetary wakes is N(e), so by comparison with 
eq. (8), 

4 Z v (e) , x 

Despite the relations (3) and (5), we are not yet in a position to evaluate a e g numerically because 
we have not determined the width of the distribution <p, or equivalently, the damping length. 

In the rest of the paper we assume a minimum mass solar nebula MMSN, cf. Hayashi 1981: 

E(r) = 1700 g cnT 3 r2u 2 , (H) 
c(r) = 1.2 km s _1 r^/ 4 , (12) 

where tau is the radius r expressed in astronomical units. For these parameters, 

Q = 67 rfj\ (13) 

which is effectively infinite for our purposes, so that we completely neglect self-gravity of the disk. 
The MMSN plausibly corresponds to the state of protostellar disks as observed around T Tauri 
stars, but disks are likely to have been more massive at earlier evolutionary phases when the star 
was enshrouded by infalling dust and gas. Self-gravity could then have been important to the 
effective viscosity and accretion rate (Gammie 2000). 

Another important assumption is that the sound speed is independent of height above the 
midplane. This is appropriate for the later evolutionary phases of protostellar disks, when the local 
thermal structure is probably dominated by passive reprocessing of radiation from the central star 
(Chiang & Goldreich 1997), and it is in these phases that planets are most likely to be dynamically 
important. Vertical isothermality allows us to treat even short-wavelength acoustic disturbances 
two-dimensionally, i.e. in radius and azimuth. When the disk is prodominantly self-luminous, as is 
likely in earlier phases and during FU Orionis outbursts, the temperature may decrease away from 
the midplane, in which case acoustic waves will concentrate toward the surface and perhaps damp 
more quickly than we estimate here (Lubow & Ogilvie 1998; Ogilvie & Lubow 1999). The thermal 
stratification of an active disk is very uncertain, however, since the distribution of dissipation over 
height is poorly understood; in some models, dissipation occurs mainly near the surface, which 
would also tend to make the disk vertically isothermal (Gammie 1996). 
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3. Linear damping 

Linear damping mechanisms are those for which the damping length is independent of am- 
plitude and which do not couple different azimuthal harmonics of the wake. One possibility is a 
background viscosity vq that acts on density waves as well as the mean shear. Takeuchi et al. (1996) 
showed that in this case the wave's angular momentum decays with distance as (for definiteness, 



consider r > r p ) fj jTn oc exp 



-2 J Ak u dr 



, where fj m is the angular momentum flux carried by 



the m th azimuthal harmonic, and Ak v is the contribution to the imaginary part of the WKBJ 
wavenumber due to viscosity: 

Mkf „ " 2 ("-y 2 -* 2 . (14 , 



c 

and 

'4 K 2 



Aki PS Vq 



— + 



=%^W). (15) 



3 m 2 (n-n p ) 2 

The most important contributions to the flux are for m ~ fi, where \i ~ r/l 3> 1 is the torque 
cutoff (GT80), so the terms involving the epicyclic frequency k (~ Q) can be neglected beyond a 
few I from the planet. Hence Aki ~ (r — r p ) 2 /l 4 , so 

fj, m -f { Ze W [-a (r-r p f/h"]. (16) 

Here = u G /VLh 2 and h = c/VLr = (3/2)/ is the disk thickness. Hence the damping length is 
Id ~ «o X ^ 3 ^- ^ a o is as large as 10~ 3 , then there is no need to invoke density waves to explain 
accretion. On the other hand, if ao 10 -3 then Id » Wh. 

Cassen & Woolum (1996) have proposed another linear damping mechanism: radiative losses 
from the disk surface. They conclude that the effect is strong, but we believe that they made 
inappropriate use of the adiabatic approximation. Cassen & Woolum estimate the radiative losses 
by evaluating the adiabatic (lossless) density-wave eigenfunction at the disk photosphere. Because 
they take an isothermal background state, as appropriate for a passive disk warmed by the central 
star, the first-order temperature perturbation at the photosphere is as large as it is at the disk 
midplane, so that their damping rate is almost independent of the total optical depth of the disk 
(cf. their Fig. 1). This result is implausible. By contrast, we find that if the non-adiabatic effects 
are included self-consistently in the linear analysis, then the damping rate is inversely proportional 
to the total optical depth and is very much reduced compared to their estimate. The typical optical 
depths in the MMSN are r ~ 10 3 - 10 4 . 

For the planetary wake, radiative diffusion in the radial direction [neglected by Cassen & 
Woolum] is much more important than losses from the surface, because the radial wavelength is 
typically shorter than the vertical thickness. The process is similar to damping of a sound wave in 
air by thermal conduction. Adapting formulae from Landau & Lifshitz (1959), we find 

k 2 h 2 2m H a B T 3 

Akl ~ - — mzr> (17) 
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where gb is the Stephan-Boltzmann constant, 2mn is the mass of the hydrogen molecule, and ks 
is Boltzmann's constant (and c is still the speed of sound, not of light). Using eq. (14) for k, we 
find that the damping length due to radiative diffusion is 

, / n^k B T \ 1/3 , , 

^HwJ • <18> 

The most interesting feature of this expression is that now lj oc r 1 / 3 , so that it is much smaller 
than in the case of radiative losses from surface where (in our analysis) Id oc r. Numerical estimates 
show, however, that the damping length is still > r for typical conditions. 

We conclude that the obvious linear damping mechanisms produce damping lengths as large 
as the radius itself, or larger. Therefore, these mechanisms are unimportant compared to shock 
damping which, as will be seen, leads to a damping length that scales with planetary mass as 

—2/5 

oc M p but is still only a modest multiple of the disk thickness even if M p is as small as M®. 



4. Linear excitation of a planetary wake 

Even if the planet does not open a gap in the disk, there is a minimum radial separation 
\r — r p \ = I at which density waves can be excited because at smaller distances the velocity of the 
gas in the corotating frame of the planet is subsonic, and a stationary perturber does not excite 
acoustic waves in a subsonic flow (Landau & Lifshitz 1959). This explains the torque cutoff (GT80), 
at least in nonselfgravitating disk where the density waves are basically acoustic. For a planet of 
sufficiently low mass [eq. (31)], the forces exerted at this minimum distance are weak enough so that 
the excitation of the wake and its angular momentum flux can be accurately calculated by linear 
theory. It is well known that in the absence of linear damping (viscosity, thermal conductivity), a 
sound wave of small but finite amplitude propagating into still air eventually shocks; the distance 
to the shock is proportional to the wavelength and inversely proportional to the initial amplitude 
(Landau & Lifshitz 1959; Whitham 1974). In a disk the shock length is shorter and the dependence 
on initial initial amplitude is weaker because differential rotation compresses the radial width of 
the wake as it propagates. 

Linear excitation has been studied by Goldreich & Tremaine (1978, henceforth GT78), GT80, 
Artymowicz (1993), Ward (1997), and others. The relative phases of the azimuthal Fourier har- 
monics of the wake are normally discarded because they are not relevant to the torque. But the 
onset of the shock depends upon the local slope of the wake front, which in turn depends upon 
these phases. Therefore, we have repeated the linear calculation using the methods of GT78 and 
GT80. 

As usual, x = r — r p and y = r p (6 — 6 p ), denote pseudo-cartesian radial and azimuthal 
coordinates in a corotating system centered on the planet, and all properties of the background 
flow are expanded to lowest order in x/r p . The background surface density £ and sound speed c 
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are made constant, as are the rotation rate (=Coriolis parameter) Q,, shear rate 2A = rdQ/dr, and 
vorticity 2B = 2(Vl + A). The background flow with respect to the planet is 2Axe y . No distinction 
is made between Q and the angular velocity of the planet, since this is not important for the angular 
momentum flux carried by the wake to leading order in h/r p [but vital to the net torque on the 
planet (Ward (1997)]. Although we are interested in keplerian disks, for which 2A/Q = —3/2 and 
2B/Q = +1/3, the physics is clarified by retaining A and B in basic formulae. We use the Mach-1 
distance I = c/\2A\ rather than h = c/VL = (3/2)1 as a reference length. A convenient unit for the 
planetary mass is 

M l = (19) 
\2A\G v ; 

which reduces to 2c 3 /30G in a keplerian disk. In linear theory, the amplitude of the wake is oc M p , 
so it is sufficient to calculate as if M p = Mi and scale the results accordingly. Note that M p > Mi 
is also the point at which the linear approximation begins to fail; the Roche lobe of the planet 
becomes > and a gap may open in the disk (§6). 

The wake is steady in the planet frame, so that its perturbations in radial and azimuthal 
velocity, u,v <C c, and the relative perturbation in surface density, a = 5S/S <C 1, are functions of 
(x,y) only. Their spatial Fourier transforms u, v, and a (k x ,k y ) satisfy (GT78, GT80) 

^ + [c 2 k(r) 2 + k 2 ]v = -ik y ^ + 2ik x (T)B0, 

ClT^ (IT 



The timelike variable r is related to the pitch angle of the wave: 



^ + 2Bk x v + ik 2 y ^j . (20) 



the notation is standard and will not be confused with optical depth. An individual Fourier compo- 
nent evolves at constant k y but varying k x due to the differential rotation. Since A < 0, all Fourier 
information flows along the characteristics dk x /dr = —2Ak y from leading (r < 0) to trailing (r > 0) 
waves, and eqs. (20) describe the evolution of the amplitudes along those characteristics. The quan- 
tity k is the instantaneous wavenumber \Jk% + k 2 , and <f> = —2-KGM p /k is the Fourier transform 
of the point-mass planetary potential. Appropriate initial conditions are i = at r = - oo. Unlike 
GT80, we have neglected the selfgravity of the disk. 

Unless otherwise noted, the rest of this section assumes units c = — 2 A = 1, so that 1 = 1, and 
the wake amplitude is scaled to M p = Mi. 

We solved eqs. (20) numerically to obtain (u, v, a) on a grid of N x x N y = 4096 x 8192 points 
in (k x ,k y ). The very large grid proved necessary to minimize numerical artifacts of the discrete 
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Fourier transform (DFT). We chose the maximum (Nyquist) value of k y = 8, well beyond the 
torque cutoff. The corresponding azimuthal grid resolution and periodicity length are Ay = (tt/8) 
and N y Ay ~ 3217, respectively. Since the locus of the wake is approximately y ~ — x 2 sign(x)/2, 
we chose L x = 2 A /L^ ps 113., implying a radial resolution Ax = L x /N x = 0.028. 




The Fourier data were filtered in pitch angle to avoid aliasing. When k x S> k y & k, the 
homogeneous solutions to eqs. (20) change phase by k x Ar = rAk x from k x to k x + Ak x , where 
Ak x = 2tt/L x is the resolution in Fourier space, and k x = rk y in the present units. This phase 
change should be < it, in other words r max = L x /2 ps 57, else the phase will be aliased and the 
inverse DFT will assign the wave to the wrong x position. We multiplied the Fourier components 
by the pitch-angle filter 



Because the pitch angle increases linearly with |x| in coordinate space, this filter attenuates the 
wake at |x| > L x /4. We also used a similar filter in the k y direction, tapering linearly from unity 

at \ky\ = to zero at \k y \ = fc^Nyquist- 

We checked our calculation against the angular momentum flux reported by GT80 in two 
different ways. One can show that as L x , L y — > oo, 



the average being taken over oscillations in r at fixed k y . We evaluated a numerical approximation 
to this formula in Fourier space (before applying the windows discussed above) and found agreement 
with the flux (3) reported by GT80 for Q = oo to the accuracy they quoted, viz. two significant 
digits. We also evaluated the flux as a function of x in coordinate space using eq. (32). The 
errors are larger because the coordinate-space wake is softened by the pitch-angle and k y filters, 
but nevertheless fj(x) is nearly constant and correct within ±3% between |x| = 6 and |x| = 28. 

Fig. 1 shows that beyond a few Mach lengths (I) from the planet, the linearized wake profile 
is approximately constant as a function of y + x 2 sign(x)/2, as expected if the source terms are 
negligible at |x| 3> 1 and the wake propagates as a tightly-wrapped acoustic disturbance (see 
below). To the extent that this is true, the profiles can be regarded equally well as azimuthal cuts 
across the wake at fixed radius, or as radial cuts at fixed azimuth but with x 2 /2 rather than x as 
radial coordinate. If plotted against x, the profile would shrink in width cx |x| _1 for \x\ S> 1. The 
profile is about as compact as possible given the torque cutoff [the integrand of eq. (21) peaks at 
k y ps 0.36 showing that the phases of the Fourier components of the wake are closely aligned. 




'1 if r < r max /2, 

W(k x , k y ) = I 2(1 - r/r max ) r max /2 < r < r, 

,0 T J> T max . 



oc 




(21) 



o 
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Fig. 1. — Surface-density profile of the linearized wake at x = —2 (heavy line), —4, —6, —8. Division 
by \x\ l l 2 removes growth due to flux conservation. Arrow shows direction of gas flow relative to 
profile. See text for dimensional units. 
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5. Nonlinear propagation 

The fully nonlinear steady-state fluid equations in xy are 



(2Ax + v)d y u + ud x u — 2VLv + —d x T, 



(2Ax + v)d y v + ud x v + 2Bu + j;d y Il 



(2 Ax + v)d y T, + ud x Y, + Y,d x u + T,d y v 



0, 

0, 
0. 



(22) 



Henceforth £ denotes the full surface density, and So its unperturbed value. The sound speed 
follows the adiabatic law c 2 (S) = Co(£/£o) 7_ \ a good approximation for weak shocks. We have 
dropped the planetary source terms, since at distances 3> I they are unimportant. 

The system (22) is hyperbolic, and in the Appendix, we demonstrate that for |x| S> I, it can 
be reduced to a single first-order nonlinear equation: 



dtX ~ xd v X = 0, 



(23) 



which is the inviscid Burger's equation (Whitham 1974). The dimensionless variables appearing 
here are related to radius, azimuth, and density contrast as follows: 
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(24) 
(25) 
(26) 



5.1. The shock 

For almost any smooth choice of initial conditions in Burger's equation (23), the solution will 
eventually threaten to become double- valued, so that one must fit in a shock (Whitham 1974). The 
identity (dtx)r]/(d v x)t = ~(9tV)x converts eq. (23) to the linear equation 

which says that each point on the wave profile moves forward at a "speed" determined by the value 
of x at that point. A shock develops when higher (larger x) parts of the profile overtake lower ones. 
The first shock develops where the initial profile is steepest, after a delay 

(28) 

o 



^shock ^0 



max 



dx ■ 

ar/ SlgnX 



-13- 



0.5 



x 








-20 



10 





7) 







20 



Fig. 2. — Nonlinear evolution of the wake profile according to eqs. (23)-(26), integrated via second- 
order advection scheme (van Leer 1977). Initial conditions taken from linear wake profile at xo = 
—21 (t = 1.903) scaled to planet mass M p = Mi. Profiles shown at t — to = 0, 4, 16, 64, 256 (highest 
to lowest). The latest profile has the expected N-wave shape. 
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If the initial profile changes sign, as it does in Fig. (1), there will be at least two shocks eventually: 
one moving to the left and one to the right in this reference frame, which is moving radially away 
from the planet at Cq. 

The shock jump condition depends upon what is to be conserved. In our case, x is approxi- 
mately proportional to the density contrast (E — So)/So provided that the latter is small, so mass 
conservation demands 



where, xi an d X2 are the values of % before and after the shock, and U s = (c^/cfa) shock is the 
"velocity" of the shock in these variables. In physical variables, the radial velocity of the shock is 
the mean of the sound velocities on either side of the jump, as usual for weak shocks. 

For an initial pulse that changes sign and has finite range, the final asymptotic behavior is 
an "N-wave," called such because the profile resembles that letter of the alphabet (Fig. 2). The 
unperturbed gas flows from left to right in these variables, so it first encounters the lefthand shock, 
followed by a rarefaction wave in which the pressure falls below its unperturbed value, and finally 
a jump back to the unperturbed conditions at the righthand shock. The point where x = ( on the 
rarefaction wave) does not move, and the flux vanishes there, so mass conservation implies that the 
areas in the left and right lobes of the N-wave are both constant. From this and the jump condition 
(29), it follows that the height and width of the N-wave scale as i -1 / 2 and i +1 / 2 , respectively (cf. 
Whitham 1974; Landau & Lifshitz 1959). 

We choose xq = —21 as the matching point between the linear and nonlinear calculations. 
This is not far enough from the planet to be completely outside the excitation region: the flux in 
the linear wake is about 7% smaller than its asymptotic value. For M p <C M\ one should probably 
use a larger value of |xo|. In practice, the interesting values of M p will be not much smaller than 
Mi, and it is then necessary to match at a rather small |xo| in order to remain within the range 
where the linear calculations are valid. 

Applying the criterion (28) to the linear wake profile at |xo| = 2 1 [see Fig. (1)], whose amplitude 
is oc M p , we find that the shock develops at t s hock — ~ 0.79(Mi/M p ). The constant to ~ 1-89 
is oc |xo//| 5 ^ 2 but independent of M p . In order that subsequent formulae have a simple power-law 
scaling with M p , we shall neglect to, but it becomes important as M p approaches Mi. Using the 
relations (24) and h = c/fi, = (2/3)/, t s hock translates to 



This result indicates that if M p > Mi, then the shock begins immediately and cannot be separated 
from the excitation region. Note that for M* = M Q , 



(Xi - U s )xi = (X2 - U s )x2 u s = 



Xl +X2 

2 



(29) 




(30) 
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5.2. Angular momentum flux and energy dissipation 

By summing over the contribution of each azimuthal harmonic as given by GT80, one can 
show that the angular-momentum flux carried by the wake in a nonselfgravitating disk is 



n/r p 

|2Ax|S 



i^rfe" / (S-S ) 2 dy. (32) 



-7r/r p 

In terms of our dimensionless variables, 

where *(«) , („,*)*,. (33) 

We have used the fact that the wake depends upon radius and azimuth most rapidly in the combi- 
nation rj, and that its slowly- varying amplitude can be thought of as a function of radius or azimuth 
according to the correspondence \x\ « \f2l\y\ = 2 1 / 2 (5t/4) 2 / 5 /. Without shocks, the dimensionless 
flux $ is strictly conserved; in fact, the integral over r\ of any function of x ls conserved. 

Once the shock develops, fj decays with distance as the angular momentum carried by the 
wake is transferred to the mean flow. From the scalings for the N-wave discussed above, it follows 
that <J>(i) oc t~ l l 2 at large t, hence 

fj(x) oc |x|- 5 / 4 (M > M shock ), (34) 

which was also confirmed by direct analytical N-wave expansion of the original set of equations 
(22). The full evolution of <J>(i) is shown in Fig. 3. 

The total energy dissipation rate associated with the absorption of the wake is [cf. eq. (4)] 

oo 

dfj 



= 2 J [n(r) - n( 



3/5 



5/ \l + 



2 00 



where we have used integration by parts to transfer the radial derivative from fj to O(r). Since 
£-3/ 5 (|>(i) oc t -11 / 10 at large t, the integral converges, but slowly; its value is ~ 12.3 when M p = M\. 

By scaling our numerical results to other M p and 7, we find that the dissipation produced by 
a small planet is 
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Fig. 3. — Solid line: dimensionless angular momentum flux [eq. (33)] carried by the wake versus 
propagation pseudotime (24), based on numerical integrations shown in Fig. 2. Dashed: Asymptotic 
t -1 / 2 scaling appropriate to the N-wave, normalized to agree with the constant preshock value 3>(io) 
at t = tshock (bullet). 
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in a keplerian disk. From eq. (10), it now follows that effective viscosity parameter is 




(36) 



For simplicity we have taken all planets to have the same mass in the final line above, and we have 
used the MMSN model (11)-(12). 

In the Abstract, we stated that shock damping is local but only weakly so. It is local in the 
sense that most of the contribution to e would come from a finite range of \x\/h in the limit of a 
very thin disk, r/h — > oo. But quantitatively, since protostellar disks are not extremely thin, the 
results of this section show that the dissipation will extend over a range comparable to the radius. 
Although fj(x) falls to half of its full value (3) at \x\ ~ 2.1 |x| s hock) half of the dissipation occurs 
beyond 9.1 |x| s hock> and the residual contribution from distances > \x\ falls off as |x| _1//4 when 

| X | S> | X | shock- 



We have shown that density-wave wakes excited by planets in a nonselfgravitating protostellar 
disk serve as an effective viscosity mechanism, provided that the waves are damped over a distance 
smaller than the orbital radius of the planet. The effective viscosity is proportional to the damping 
length, which is most likely to be limited by shocks unless the mass per planet M p <C M®. If a 

3/5 /i 

fixed total mass is available to form planets, then a e g oc M p . Eq. (36) suggests that a c g ~ 10 
in the inner solar system where M p ~ M®, and a e s ~ 3 x 10 -4 at r ~ 10 AU if M p ~ 10 M®, 
which is characteristic of the rock+ice cores of the giant planets, though very much less than the 
total mass of Jupiter or Saturn (Guillot 1999). 

Our estimates for a e g may seem optimistic because they assume that all of the rocks and 
ices are efficiently converted into planets. On the other hand, a c g could be even larger if M p is 
augmented by gas. In that case, a e s oc M p ^ 3 because the mass in planets is not limited to the 
metals. There are two restrictions on the amount of gas we can allow for. The first is technical: 
we have calculated the excitation of the density waves from linear theory, which cannot be justified 
when M p exceeds the limit (31). This could be overcome by nonlinear two-dimensional numerical 
calculations (e.g. Korycansky & Papaloizou 1996). The other restriction is more substantial: once 
M p is large enough to open a gap in the disk, the dissipation rate probably scales more slowly than 



Several estimates exist for the minimum mass required to open a gap, M gap . Lin & Papaloizou 
(1993) and Bryden et al. (1999) argue that in a completely inviscid disk, M gap « M\ as given by 
eq. (31) because of radial gas pressure gradients; recall that at this mass, the Roche lobe (or Hill 



6. 



Discussion 
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sphere) of the planet is h. Ward (1997) argues for M gap » 0.147r£/ir, or M gap » 0. 8 (r /At/) 3/4 M® 
in the fiducial MMSN (11)-(12). This results from comparing the time to open a gap with the time 
to migrate across it, assuming that planetary wake damps immediately at the resonances where it 
is excited. Presumably M gap should be increased to allow for the finite damping length that we 
have computed, but it would still be proportional to S, unlike eq. (31). 

Mg a p is expected to be larger in a viscous disk than in an inviscid one because the torque 
exerted by the planet can be balanced by a viscous torque. Bryden et al. (1999) cite 



In our case, it seems likely that the shocks created by one planet will resist the opening of a gap 
next to an adjacent planet, so it is perhaps reasonable to use the value (36) for a on the righthand 
side above. 

In the discussion to follow, we shall follow Lin & Papaloizou (1993) and adopt (31) as a 
provisional value for M gap , but clearly the question deserves further study. Then M gap ~ 45 M® 
at 10 AU using eq. (12). If the planets at 10 AU accrete enough gas to reach this limit (Saturn 
exceeds it by a factor ~ 2), then a e g ~ 3 x 10~ 3 . 

In short, values for a e s in the range from 10 -4 up to at least 10~ 3 seem plausible for the 
density-wave mechanism in a minimum-mass solar nebula. The predicted accretion timescale is 



in reasonable agreement with the disk lifetimes cited in §1. 

By comparing observed sizes, ages, and accretion rates of T Tauri disks with self-similar models, 
Hartmann et al. (1998) estimate a ~ 10 -2 , which may be too large to be explained by planetary 
wakes, especially if one must rely on the rocky core masses only. The strongest lower bounds on a 
assume that the disks started out much more compact than observed. In view of the wide scatter 
and uncertainties in the observationally derived parameters, the questionable relevance of self- 
similar models, and the theoretical uncertainties just discussed, it is unclear whether Hartmann 
et al.'s result is a major difficulty for us. Another concern is that planets are likely to migrate 
inwards with respect to the gas (Ward 1997), perhaps enhancing the viscosity of the inner disk but 
diminishing that of the outer parts where most of the gas probably resides. Still another is that our 
mechanism operates only after planet growth, which standard models place late in the lifetime of 
the gaseous nebula, if not afterwards (Podolak et al. 1993; Ikoma et al. 2000), whereas Hartmann 
et al. find larger accretion rates in younger disks on average. 




which is again independent of S (if a is). In the MMSN, this becomes 
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All of these objections have merit. Nevertheless, the density-wave mechanism is interesting 
from several points of view. First, it provides a minimum effective viscosity when other mecha- 
nisms are in abeyance. For example, MHD turbulence may dominate only during FU Orionis-like 
outbursts when the temperature and ionization of the disk rises (Gammie 1998); self-gravity may 
dominate in young, massive, cool disks (Lin & Pringle 1987); but neither is likely to be effective in 
the late, thermally passive phases when it is believed that planets form. Second, calculable models 
for the effective viscosity are so difficult to come by that it is worthwhile to examine any plausible 
candidate, even if others seem more likely to dominate (but cannot be calculated). Third, the 
formalism offered in §5 and the Appendix for the development and damping of the shock, though 
only approximate, is expected to be reasonably accurate for weak shocks and simple enough to 
use that it may find other applications. 1 Finally, there may be other nonaxisymmetric structures 
in protostellar disks of greater mass than planets, for example vortices (Adams & Watkins 1995; 
Godon & Livio 2000; Nauta 2000). Any such structure would create a density- wave wake. 

We thank Charles Gammie, Peter Goldreich, Brad Hansen, Kristen Menou, and Scott Tremaine 
for helpful discussions. This work was supported by the NASA Origins Program under grant NAG5- 
8385. 
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A. Reduction to Burger's equation 

At \x\ I, eqs. (22) are hyperbolic, as usual for supersonic steady flows (Landau & Lifshitz 
1959): that is, they describe an initial- value problem. It is convenient to regard y as the timelike 
variable. The three characteristic "velocities" are 



dx\ _ u u(2Ax + v)± c y / (2Ax + v ) 2 + u 2 - 



2 



- 



dy) + _ 2Ax + v } (2Ax + v) 2 -c 2 

Eqs. (22) simplify for \x\ » /, where a linearized WKB treatment suggests that d y <C d x and d«u. 
Dropping v entirely, and d y except in the combination 2Axd y , results in 

c 2 (T,) 

dyU + Ud ( U+-^-d^ = 

d y T, + ud^T, + Ed^u = 0. (Al) 

We have replaced x with £ = Ax 2 so that the equations are autonomous, i.e., the independent 
variables do not appear explicitly. 

The system (Al) is formally identical to one-dimensional isentropic gas dynamics and can be 
solved by the same methods (Landau & Lifshitz 1959). It has characteristics C± and Riemann 
invariants R±, 

C±:(^f) =±c+u, R±= u ±-^-. (A2) 

\dyj± ' 7-1 

The C+ characteristics propagate toward the planet from the unperturbed disk, where u = and 
c = Co- Therefore R + = 2co/(7 — 1) everywhere, and u = 2(co — c)/(7 + l)- 2 Eqs. (Al) then reduce 
to a single first-order equation, 

o = d y 4> - (i + V)co%v, ( A3 ) 

if c\ - 7 + 1 c - cp 7 + 1 g - S 

^(y,0 = - — 7—— y — ifV<i, 

7 — 1 Co 2 Z^o 

which is the well-known inviscid Burger's equation: the simplest nonlinear equation capable of 
displaying a shock (Whitham 1974). 

Unfortunately, the analysis above is oversimple. In the limit of an infinitesimal disturbance, 
eqs. (Al) predict that the wave propagates at constant amplitude. But a linear WKB treatment 
of eqs. (22) shows that the amplitude of such a wave grows in proportion to | a; | 1 / 2 oc as a 

consequence of conservation of angular momentum flux. To recover this behavior, the larger sub- 
dominant terms involving v that were omitted in passing from eqs. (22) to (Al) must be reinstated. 
The equations for the C± characteristics become 

d±R± ~ (2fiv T cdyV) = S± ' 



2 We take x < and y > 0, which describes the wake interior to the orbit of the planet. The discussion that follows 
is equally valid for x > and y < if one interchanges the roles of C+ and C_ 
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where R± have the same meaning as before, and we have introduced the abbreviations d± = 
dy + (u± c)d^ for the derivatives along C±. The source terms S± are smaller than individual terms 
in R± by 0(£ _1 ) but nonzero, so that the R± are not invariant. Yet R + is still well approximated 
by its undisturbed value 2co/(7 ~~ 1) : First of all R + is truly undisturbed until the incoming C + 
characteristic encounters the leading edge of the wake, which propagates in the opposite direction 
on a bundle of C_ characteristics of roughly constant width A£. Then, since the encounter lasts a 
"time" Ay ~ A£/2c, the maximum change in R + relative to its undisturbed value is 0(A£/£) <C 1, 
which may be neglected at large £. 

The source term 5_ cannot be neglected because it "travels with" the wake so that its effects 
accumulate. We solve for S- using the second of eqs. (22), which, to an adequate approximation, 
reads 

a ' v K -ib ( 2Bu + ■ 

Only those parts of S- that are in phase with u or S — So are important for the growth of the 
amplitude; the out-of-phase parts make a slight but negligible change to the characteristic velocity. 
The operator d y changes the phase by 90°, so from the equation above, 

1 1 / c 2 - c 2 \ 

r- [cdyv + 2Qv] iT , nhaKP = — 2Bcu + 2Q & 

2Ax 1 y Jm phase 2A£ V 7-1 / 

We approximate c by cq except in the difference c — cq and invoke the approximate constancy of 
R + to write u w 2(co — c)/(j — 1). Finally we write Y> = (7+ l)(c — co)/(j — l)co as before, and so 
obtain the corrected C_ equation in the form 

^-(l + V)co^ = -g^ (A4) 

As expected, the general solution of the linearized form of eq. (A4) is = |^| 1//4 /(C + coy), 
where / is an arbitrary function. The full equation incorporates all three processes that dominate 
propagation of the wake: 

(i) increasing radial wavenumber with distance from the planet (d x = 2Axd^)] 

(ii) increasing amplitude due to conservation of angular-momentum flux; 
(hi) nonlinear steepening. 

Fortunately, eq. (A4) can be transformed back into Burger's equation by changing the inde- 
pendent and dependent variables. First let £ — > cqI(t) — y/l) (so 77 is dimensionless) : 

Wytp - ipdvip = -j-ip. 

4y 

Since the wake is confined to \rj\ y/l, we have approximated y — rjl by y on the righthand side. 
Next, let ip = (y/iy^x to absorb the source term: 

id y x-(y/i) 1/ \d v x = o. 

Finally, introduce a new dimensionless "time" (or azimuth) variable t as in eq. (24), so that the 
dimensionless equation of motion becomes (23). 



